PROGRAM MAPS
USE POLYMORPHIC_COMPLEXTAYLOR
TYPE(DAMAP) Id, MAP
REAL(DP) L1,KICK1,L2,KICK2,LD 
REAL(DP) :: PREC = 1.D-10
type(vecfield) F
type(pbfield) H
integer mf,me
mf=20
open(unit=mf,file='results.txt')
CALL INIT(NO1=4,ND1=2,NP1=0,NDPT1 =0)     !   <------------------ init for maps in ND1 degrees of freedom


call alloc(Id,Map)
call alloc(F); call alloc(h);

KICK2=-0.9d0; L2=0.5D0

F=F.CUT.C_%NO   !  CUT THE HIGHEST ORDER TO MAKE IT CONSISTANT WITH THE POISSON BRACKET OPERATOR

call VECFIELD_QUAD(L2,KICK2,F)

Write(mf,*) "  "
Write(mf,*) " VECTOR FIELD OF A QUADRUPOLE TRUNCATED TO ORDER NO-1"
CALL PRINT(F,MF)


call PBFIELD_QUAD(L2,KICK2,H)


Write(mf,*) "  "
Write(mf,*) " POISSON BRACKET FIELD OF THE SAME QUADRUPOLE "

CALL PRINT(H,MF)

Id=1;

MAP=H*ID

Write(mf,*) "  "
Write(mf,*) " POISSON BRACKET FIELD ACTING ON IDENTITY MAP  "
Write(mf,*) "  SHOULD AGREE WITH VECTOR FIELD "

CALL PRINT(MAP,MF)


! One can get the vector field using the (=) sign or vice versa

F=h       ! h=F is also possible

Write(mf,*) "  "
Write(mf,*) " VECTOR FIELD FROM POISSON BRACKET  "

CALL PRINT(F,MF)

call kill(Id,Map)
call kill(F); call kill(h);
close(mf);
END PROGRAM MAPS

subroutine VECFIELD_QUAD(L,K,F)
USE POLYMORPHIC_COMPLEXTAYLOR
implicit none
TYPE(VECFIELD) F
REAL(DP) L,K

! DRIFT PART
F%V(1)=L*(1.D0.MONO.'01')  /SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))
F%V(3)=L*(1.D0.MONO.'0001')/SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))

! QUADRUPOLE PART 

F%V(2)=-L*K*(1.D0.MONO.'1')
F%V(4)= L*K*(1.D0.MONO.'001')


end subroutine VECFIELD_QUAD

subroutine PBFIELD_QUAD(L,K,H)
USE POLYMORPHIC_COMPLEXTAYLOR
implicit none
TYPE(PBFIELD) H
REAL(DP) L,K

! DRIFT PART
H=L*SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))

! QUADRUPOLE PART 

H%H=H%H-(L*K/2.D0)*((1.D0.MONO.'2')-(1.D0.MONO.'002'))

end subroutine PBFIELD_QUAD

